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ABSTRACT 

We present self-consistent models of gas in optically-thick dusty disks and 
calculate its thermal, density and chemical structure. The models focus on an 
accurate treatment of the upper layers where line emission originates, and at 
radii > 0.7 AU. Although our models are applicable to stars of any mass, we 
present here only results around ~ 1M Q stars where we have varied dust proper- 
ties, X-ray luminosities and UV luminosities. We separately treat gas and dust 
thermal balance, and calculate line luminosities at infrared and sub-millimeter 
wavelengths from all transitions originating in the predominantly neutral gas that 
lies below the very tenuous and completely ionized surface of the disk. We find 
that the [Aril] 7/im, [Nell] 12.8/im, [Fel] 24/xm, [SI] 25/im, [Fell] 26/mi, [Sill] 35 
/im, [01] 63/im and pure rotational lines of H2 and CO can be quite strong and 
are good indicators of the presence and distribution of gas in disks. Water is an 
important coolant in the disk and many water emission lines can be moderately 
strong. Current and future observational facilities such as the Spitzer Space Tele- 
scope, Herschel Observatory and SOFIA are capable of detecting gas emission 
from young disks. We apply our models to the disk around the nearby young 
star, TW Hya, and find good agreement between our model calculations and 
observations. We also predict strong emission lines from the TW Hya disk that 
are likely to be detected by future facilities. A comparison of CO observations 
with our models suggests that the gas disk around TW Hya may be truncated to 
~ 120 AU, compared to its dust disk of 174 AU. We speculate that photoevap- 
oration due to the strong stellar FUV field from TW Hya is responsible for the 
gas disk truncation. 

Subject headings: infrared:ISM — line: format ion — planetary systems: proto- 
planetary disks — radiative transfer 
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1. Introduction 

Circumstellar disks have been observed around stars with a wide range of masses, from 
brown dwarfs to massive B stars (e.g., Schreyer et al. 2006, Luhman et al. 2007). Disks 
facilitate accretion and angular momentum transport during the star formation process and 
provide a mass reservoir of gas and dust for the formation of planetary systems. Planet for- 
mation in disks appears to be fairly common, with observed extrasolar planets now totalling 
over 250 objects (e.g., Udry et al. 2007). The bulk of the initial gas and dust mass in disks 
around solar mass stars has a finite lifetime of order ~ 1 — 10 Myr (e.g., Haisch et al. 2001, 
Sicilia-Aguilar et al. 2006), many orders of magnitude smaller than stellar lifetimes, and this 
sets an upper limit on the formation time of giant planets. 

Disk studies are largely based on measurements of the dust emission. In primordial disks, 
although the dust mass is much smaller than the gas mass, the dust is a strong emitter (see 
review by Watson et al. 2007). Disks transition from a younger optically thick dust stage 
characterised by strong infrared continuum excesses, to an optically thin dust stage with 
weak infrared emission (e.g., Najita et al. 2007). The census of optically thick young disks 
and disks with scant infrared excesses like the weak-lined T Tauri stars (WTTS) indicates 
that the initial optically thick epoch lasts ~ 3 Myr (e.g., Haisch et al. 2001), the transition 
phase is rather abrupt (< O.lMyr, see Hillenbrand 2007), and that the WTTS phase also 
lasts for ~ a few Myr. The missing small dust particles in WTTS and/or transition disks may 
have coalesced to build planetary objects (e.g., Quillen et al 2004). However, observational 
studies of disks infer initial masses (~ 0.1M Q or 100 Mj) much greater than that of the mass 
of planets around stars (< 10 Mj), suggesting that much of the gas and dust is removed 
from the disk either by accretion onto the central star or by disk dispersal processes. 

The limited observations of gas in disks seem to broadly agree with the evolutionary 
sequence traced by dust emission, but dispersal timescales are less constrained at < 10 Myr 
(e.g., Zuckerman et al. 1995, Hollenbach et al. 2005, Pascucci et al. 2006). Gas giant planets 
must necessarily form on these timescales before the gas disk disperses. As gas dominates 
the mass of giant planets and determines the dynamics of solid particles (until gas mass 
becomes smaller than particulate mass), studies of the disk gas distribution and evolution 
are critical to understanding both planet formation and disk destruction processes. 

Gas disk observations are difficult because of the inherently weak line emission of the 
primary constituent, H 2 . As a result, low abundance tracers are generally used to infer the 
presence of gas. Gas emission lines have been observed from numerous nearby disks, mainly 
millimeter and sub-millimeter wavelength rotational transitions of 12 CO, 13 CO, HCO + , DCO, 
HCN, CN, near-infrared vibrational lines of H 2 and CO, and more recently mid-infrared H 2 
pure rotational lines, as well as [Fel] and [Nell] emission (e.g., review by Dutrey et al. 2007, 
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Pascucci et al. 2007, Lahuis et al. 2007). The excitation of these tracers is sensitive to 
the density and temperature of the gas as well as the chemical abundance of the emitting 
species. Therefore, diagnosing gas disk parameters from the measured line intensities requires 
detailed models of the thermal, density and chemical structure of gas in the disk. 

Detailed models of gas disks are recently becoming available in varying degrees of so- 
phistication (e.g., Kamp & Dullemond 2004, Gorti & Hollenbach 2004, hereafter GH04, 
Jonkheid et al. 2004, Glassgold et al. 2004, 2007, Nomura & Millar 2005, Aikawa & Nomura 
2006, Nomura et al. 2007). Predicted line emission from these models can be compared with 
observations to infer the physical conditions in disks. Many existing gas disk models solve 
for the chemistry and gas temperature, but assume a disk vertical density structure that is 
derived from modeling the dust radiative transfer (e.g., Glassgold et al. 2004, 2007, Semenov 
et al. 2006, Meijerink et al. 2007). Kamp & Dullemond (2004) calculated the gas density and 
temperature by adopting the disk structure determined by the dust models of Dullemond et 
al. (2002), whereas Jonkheid et al. (2004) and Glassgold et al. (2004, 2007) use the dust 
disk model of D' Alessio et al. (1998). Because the gas temperature may differ from the dust 
temperature and the gas/dust ratio may change with vertical height, these calculations may 
be inconsistent with the condition of vertical pressure equilibrium. These inconsistencies are 
especially severe in the upper layers of disks where gas and dust temperatures can deviate 
significantly. Moreover, most observed gas emission originates from these layers making any 
interpretation of data contingent on obtaining the correct gas temperatures and densities at 
these heights. GH04 and Nomura & Millar (2005) address this issue and self-consistently 
calculate the gas density and temperature as a function of disk radius and vertical height. 
In these models, the dust temperature is separately determined and the gas temperature 
obtained from a detailed thermal balance and chemistry calculation, by considering various 
heating and cooling sources. GH04 only consider optically thin dust disks, but solve for the 
gas disk in detail with heating by dust collisions, X-rays, grain photoelectric heating due to 
far ultraviolet (FUV) radiation incident on very small grains, cosmic rays and exothermic 
photoreactions, and cooling by many atomic, ionic and molecular transitions. Nomura & 
Millar (2005) consider younger disks with optically thick dust, but their gas model is more 
simplistic. They consider heating by dust collisions, cosmic rays, and grain photoelectric 
heating, and cooling is only by OI, CII and CO lines. Aikawa & Nomura (2006) recently 
extend these models to compute detailed disk chemistry, but use a disk structure (density 
and temperature) derived from the earlier simpler models, and Nomura et al. (2007) include 
the effects of X-ray heating. 

Our earlier models (GH04) of gas emission from evolved, optically thin dust disks com- 
bined with recent Spitzer observations of disks around stars of various ages (FEPS Legacy 
Science Program), found that disks with optically thin dust also have very little gas (Pascucci 
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et al. 2006). We could set stringent limits on the gas mass using our models in conjunction 
with this data and estimated upper limits on gas disk lifetimes to be 5 — 30 Myrs. This result 
suggests that as inner (~ 1 — 20 AU) disks transition from optically thick to optically thin, 
they also lose much of their mass in gas. The same result was independently determined by 
a recent study of a large sample of disks in Chameleon, where disk gas accretion rates were 
found to be correlated with the presence of an inner dust disk, suggesting that dust and 
gas are removed simultaneously from the inner disk (Damjanov et al. 2007). Gas dispersal 
therefore presumably begins at earlier stages of evolution when the disks are still optically 
thick in dust. 

In this paper, we extend our earlier detailed gas disk models (GH04) to consider younger 
disks with optically thick dust. We are primarily concerned with calculating line emission 
from young disks that can be used as gas diagnostics to trace its distribution and structure 
at different evolutionary stages during its optically thick dust epoch. We also propose to use 
these disk models to compute the mass loss rate due to photoevaporation from stellar EUV, 
FUV and X-ray photons in an accompanying paper. Our models are constructed with the 
necessary physics to determine the disk thermal and chemical structure in the surface layers 
of the disk where line emission arises and photoevaporative flows originate. 

The structure of the paper is as follows. We describe the details of our disk models in 
§2. We discuss the resulting disk structure and line emission in §3. In §4, we apply our 
models to observational data from the well-studied disk of TW Hya and discuss the results. 
We summarize and conclude in 55. 



2. Model description, inputs and assumptions 

2.1. Main model features 

The disk models presented here are an extension of GH04. Our earlier models were 
restricted to cases where the dust was optically thin to stellar radiation, keeping our dust 
radiative transfer simple. However, the gas transitions in these models could be optically 
thick or thin depending on gas column density, and gas radiative transfer was calculated 
using an escape probability formalism. In this paper, we generalize our models to consider 
optically thick dust and therefore allow modeling of younger disks which are opaque in the 
midplane regions to stellar visible photons. In most of our models, disks are also opaque to 
their own infrared photons. Other improvements are the addition of new coolants such as 
fine structure lines of Ne and Ar ions and vibrational lines of CO, heating by and chemistry 
of Polycyclic Aromatic Hydrocarbons (PAHs), a mixture of dust chemical compositions, and 
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a model for the EUV heated surface layer. We treat in a separate paper the emission of fine 
structure lines such as [NeII]12.8/im from the EUV-heated and completely ionized surface 
layer, and here only use the model to determine the vertical position of the ionization front 
that separates this layer from the predominantly neutral disk underneath it (Hollenbach 
et al. 1994). In this paper, we do, however, treat the ionized fine structure lines such as 
[NeII]12.8/xm produced in the "neutral disk" where X-rays partially ionize and heat the gas. 

We briefly summarize the disk model here and in the remainder of this section provide 
more details for the interested reader. Our disk models solve for chemistry and thermal 
balance, and impose vertical hydrostatic equilibrium to separately calculate the density 
and gas and dust temperatures as a function of spatial location in the disk. We consider 
heating of the gas due to X-rays, grain photoelectric heating by PAHs and small grains, 
cosmic rays, exothermic chemical reactions, formation heating of H 2 , collisional de-excitation 
of vibrationally excited (by FUV) H2, photodissociation and photoionization of molecules 
and atoms, and collisions with warmer dust grains. Cooling of gas is by line emission 
from atoms, ions and molecules, and by collisions with cooler dust grains. Our chemical 
network is moderate and focused towards including species that are dominant coolants in 
the disk (enabling an accurate determination of the gas temperature). We treat gas and 
dust independently, and allow for different spatial distributions (i.e., a position-dependent 
gas/dust mass ratio), although we generally consider them well-mixed and the ratio constant 
throughout the disk. We consider a mixture of chemical compositions for the dust, and a 
range of grain size distributions. Dust radiative transfer for optically thick dust disks is 
simplified, so as to keep the numerical disk model computations tractable. We use the two- 
layer model for dust with some modifications as first described by Chiang & Goldreich (1997). 
We include the effects of background infrared radiation due to dust on line transitions of the 
gas species (see Hollenbach et al. 1991). We emphasize that our models are hydrostatic, and 
assume steady-state chemistry and thermal balance. We also neglect freezing of gas species 
on to cold dust grains. This process is likely to be important only in the dense midplane 
of the outer disk (e.g., Aikawa et al. 1999, 2002), a region unimportant for many of the 
considerations of this paper. Accretion heating can be important in the dense midplane, 
especially of the inner disk (e.g., D'Alessio et al. 1998). We ignore this term, however, 
because we are interested in evolutionary stages of the disk when the accretion rates are 
low and because gas emission lines mostly originate from the upper layers of the disk. In 
the modeling procedure, the input parameters include stellar characteristics, the initial gas 
phase abundances of the 10 elements, the surface density distribution, radial extent and 
dust properties. We then calculate the gas temperature, density and chemical structure as 
a function of spatial location (r, z) throughout the disk. 
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2.2. Stellar Characteristics 

The properties of the central star are basic inputs to our disk models. When modeling 
observations of specific star-disk systems such as TW Hya in §4, stellar properties are taken 
from the available literature. In §3 of this paper, we aim to be more general and use 
typical stellar characteristics to construct standard or fiducial disk models that depict our 
calculations of the structure and line emission. 

The models presented in this paper assume steady-state thermal and chemical balance at 
a time when the disk surface density distribution follows a prescribed power law with radius. 
They are not time-dependent calculations. However, stellar properties evolve with time, the 
disk mass and surface density distribution evolve, and the accretion rate onto the central star 
decreases (Hartmann et al. 1998) with time. Younger stars have stronger radiation fields 
which would result in greater heating and enhanced emission line strengths. For our fiducial 
disk model, we use a median set of stellar characteristics during this evolution, typical of a 
1-2 Myr old system. We use the pre-main-sequence evolutionary tracks of Siess et al. (2000) 
to establish the basic properties of the star, such as its mass, radius, bolometric luminosity 
and effective temperature. In this paper, we focus on disks around 1M Q stars. 

The radiation field created by the central star greatly affects the disk structure. Disk 
models to date have primarily focused on the dust component which absorbs the optical 
photons from the star, and optical spectra of young stars are very well determined. Gas is 
then assumed to be collisionally coupled with the dust and follow the same temperature dis- 
tribution. However, in our models, we independently calculate the gas temperature. Apart 
from indirect heating by stellar optical photons via collisions with the dust warmed by these 
photons, gas can also be heated by the ultraviolet and X-ray flux from the star. These pho- 
tons also influence gas chemistry which can affect disk structure indirectly by influencing the 
cooling which sets the gas temperature structure (see GH04 for details on these processes). 
Hence we need to characterise the stellar spectrum at these shorter wavelengths, which are 
less well determined than the optical spectra. 

X-rays Recent observational studies using ROSAT and Chandra have significantly fur- 
thered our understanding of the nature of the X-ray emission and pre-main-sequence mag- 
netic activity of young stars (e.g., Feigelson & Montmerle 1999, the COUP Project, Feigelson 
et al. 2005, Preibisch et al. 2005). The X-ray luminosity of low mass young stars is believed 
to mainly originate in the chromosphere and observed to be correlated with the bolometric 
luminosity, log Lx / 'L^oi ~ —3.6, although with a large scatter (Preibisch et al. 2005). For 
more massive stars, (M* > 3M ) the X-ray luminosity decreases and logLx / 'L^oi ~ —6 
(Preibisch et al. 2005, Steltzer et al. 2005). We therefore scale the total X-ray luminosity of 
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a star with its bolometric luminosity according to these empirical relations. X-rays originat- 
ing from accretion processes peak at lower energies, and new Chandra and XMM-Newton 
X-ray spectroscopy of a few nearby, young accreting stars reveal both a low and high energy 
component (Telleschi et al. 2006). Since the X-ray absorption cross-section increases with 
decreasing photon energy (Wilms et al. 2000), a softer X-ray spectrum would heat the lower 
density gas at the disk surface, whereas harder X-ray photons will penetrate deeper into the 
disk to heat higher density gas. We therefore extend our X-ray spectrum from 0.1 to 10 
keV. Our adopted X-ray spectrum (Giidel et al. 2007, Feigelson & Montmerle 1999, also see 
GH04) peaks at 2 keV and is approximated by a simple power-law with dLx/dE ~ E for 
0.1 < E < 2 keV and dL x /dE ~ E~ 2 for 2 < E < 10 keV. 

UV FUV (6eV< hv <13.6eV) spectra of young stars are more difficult to determine obser- 
vationally. The spectra are also highly variable, heavily extincted and of somewhat uncertain 
origin (e.g., IUE, Valenti et al. 2003; FUSE spectra, e.g., Bergin et al. 2003, Herzceg et al. 
2004). The FUV flux appears correlated with the accretion rate, suggesting that accretion 
hotspots are the main cause of the UV excess, whereas significant FUV fluxes from low accre- 
tors also indicate a chromospheric component (Calvet & Gullbring 1998, Valenti et al. 2003). 
For the purpose of this study, we follow the results of the various studies by Gullbring et 
al. (1998, 2000) and assume that the FUV originates from accretion hotspots on the surface 
of the star. We use the empirical results of Muzzerolle et al. (2003) which suggest that the 
mass accretion rate M acc for a given stellar mass scales with M%. For a 1M star of age 1 ~ 2 
Myr, we take a typical M acc ~3x 10~ 8 M Q yr _1 . We then calculate the accretion luminosity 
[L acc = GMi:M acc /2R*). The FUV spectrum is assumed to be a blackbody at a "hotspot" 
temperature of 9000 K (e.g., Gullbring et al. 1998, 2000). We independently confirmed the 
validity of the above approach by estimating the median FUV flux from stars with good 
quality IUE spectra (Valenti et al. 2003, IUE spectral atlas). In both cases, the median 
Lpuv/Lboi ratio for solar-mass stars was found to ~ 10~ 2 to 10~ 3 . We, however, adopt the 
semi-empirical approach of an accretion-rate determined FUV flux and spectrum, because 
of the uncertainties (extinction, distances, stellar mass, age) in converting the observed IUE 
spectra to a mass-dependent FUV flux. 

The EUV luminosity of young solar mass stars is very poorly constrained by observa- 
tions, but believed to be generated by the active chromospheres from these stars (e.g., Bouret 
& Catala 1998, Alexander et al. 2005). We therefore assume an EUV luminosity similar 
in magnitude to that of the chromospheric FUV and X-ray luminosities, ~ 10~ 3 L&oj. The 
photon luminosity </>euv is thus ~ 4 x 10 41 s -1 for a 1M Q star, which then determines the 
location of the ionization front at the disk surface (e.g., Hollenbach et al. 1994). 
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Stellar parameters of our fiducial disk model are listed in Table [TJ 

2.3. Dust properties 

Thermal energy exchange by collisions of gas with dust is one of the key processes in 
disks, which even at their surfaces (Ay ~ 1 to the central star) have high gas densities 
(n ~ 10 7 cm -3 ) compared to interstellar clouds. Dust grains at the disk surface absorb 
stellar optical radiation, warming this surface layer. This warm thin layer with very little 
mass then heats the disk interior through re-radiation (Chiang & Goldreich 1997, D'Alessio 
et al. 1998, Dullemond, Dominik & Natta 2001). Collisions with dust can either heat or cool 
the gas depending on whether the dust grains are hotter than or cooler than the gas. If they 
are well-coupled, gas and dust temperatures are nearly identical. However, the collisional 
coupling may be reduced by dust settling to the midplane, a process that can significantly 
decrease the local ratio of dust to gas mass (<C 0.01— the interstellar and disk primordial 
value) and the dust cross sectional area per H nucleus, <th, near the surface of the disk. In 
addition, grains may have coagulated, a process that also decreases cth- The decrease in 
collisional coupling enhances the difference between the gas and dust temperatures (GH04, 
Kamp & Dullemond 2004, Jonkheid et al. 2004) because of other possible processes that 
can affect the thermal balance of the gas. This is especially true at the disk surface where 
densities and opacities are lower. Gas heating/cooling and disk structure thus depend on 
dust properties and are affected by both grain growth and settling. 

We consider dust grains to have a range of sizes, distributed according to the MRN 
power-law characteristic of dust in the interstellar medium (n(a) oc a~ 3 ' 5 ,a m j„ < a < a max , 
where a is the grain radius). We assume spherical grains and an interstellar dust composition. 
Grain growth up to at least millimeters in size is inferred to occur in young circumstellar 
disks (e.g., van Boekel et al. 2005, Shuping et al. 2006, Muzerolle et al. 2006, Kessler-Silaci 
et al. 2006). Dust settling is also inferred with settling timescales in disks calculated to be 
rather short (~ 10 4 years at 1 AU, Dullemond & Dominik 2004), but observations indicate 
the presence of small dust even in the surfaces of evolved disks (Hernandez et al. 2005, 
Furlan et al. 2006), suggesting the operation of other mechanisms that cause shattering of 
grains and the mixing of grains from the midplane to the surface (e.g., Dullemond & Dominik 
2005). We keep the minimum and maximum grain sizes constant in this work at 50A and 
20/im respectively, which lowers the opacity (<7h) by a factor of ~ 10 from interstellar values 
(where ctr ~ 2 x 10 _21 cm 2 ). Grain growth in the disk will lead to an increase in a max in 
general. One would also expect that since small grains coagulate to form larger objects, a m i n 
might also increase. However, as pointed out earlier, observational studies infer the presence 
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of small dust grains (and PAHs) in evolved disks(e.g., Hernandez et al. 2005, Furlan et al. 
2006, Geers et al. 2006) presumably formed by a collisional cascade process which replenishes 
the small dust population in the disk (e.g., Weidenschilling & Cuzzi 1993). We therefore fix 
a m i n in our models. The important physical parameter in the models is the dust area per 
H atom, cthj which can be reduced either by increasing the dust grain sizes or by reducing 
the dust /gas mass ratio (77) in the disk. In our modeling of dust properties, we have already 
assumed some grain growth by increasing a max from the interstellar value of ~ 0.2/jm to a 
value of 20/jHi. However, we model further possible variations in cr H due to either further 
grain growth and/or settling by varying the dust/gas ratio, 77. This allows us to include dust 
settling by setting 77 to be a function of position (r, z) in the disk. However, in this paper, 
we keep the ratio 77 a constant for any given model. To summarize, we model reductions 
in opacity due to grain growth or settling processes first by increasing a max to 20/im and 
then subsequently by simply reducing the dust/gas mass ratio in the disk from its fiducial 
(and interstellar) value of 0.01. We consider one additional case with ctr ~ 2 x 10 _24 cm 2 , 
or a reduction factor of 100 from our fiducial case. Note that although modeling grain 
growth/settling by reducing 77 does not affect gas heating, temperature or spectra, which 
depend only on the gas-grain coupling (cth), the dust continuum emission would be affected 
by this procedure. When modeling specific star-disk systems, such as TW Hya (§4), we use 
a dust grain size distribution constrained by the continuum dust emission. 

Polycyclic Aromatic Hydrocarbons or PAHs are important gas heating agents in the 
FUV-irradiated ISM. They heat the gas via FUV photoelectric ejection of energetic electrons 
into the gas. PAHs are abundant in the interstellar medium and hence are likely to be present 
in primordial disks in the very earliest stages of disk evolution. The presence of PAHs in 
disk surfaces is indicated by their observation in disks around both low-mass T Tauri stars 
and intermediate-mass HAeBes (Habart et al. 2004, Acke and van Ancker 2004, van Boekel 
et al. 2004, Geers et al. 2006). It is unclear how the PAH abundance is affected by the dust 
evolution process and in view of these recent studies we assume that disks do contain PAHs, 
but in reduced abundances. We calculate the reduction in <th for our assumed dust grain size 
distribution as compared to the ISM, and scale the PAH abundance down from that in the 
ISM (8.4 x 10~ 7 , Li & Draine 2001) with the same factor. PAHs can dominate gas heating, 
especially in the presence of strong UV fields that are likely to be present during early disk 
evolution epochs. We also consider collisional heating of gas by warm PAH molecules. We 
follow Li & Draine (2001) and Draine & Li (2001) for calculating the PAH temperature and 
absorption coefficients and Weingartner & Draine (2001) for the heating due to photoelectric 
effect by small grains and PAHs. 



-10- 



2.4. Disk properties 

For our model calculations, the disk mass, gas/dust ratio, radial extent and surface 
density distribution are input parameters. Model disk masses, motivated by the above 
studies, are assumed to scale with the mass of the central star and we adopt M^k/M* ~ 0.03 
(e.g., Beckwith et al. 1990, Natta et al. 2001, Andrews & Willams 2005). We also assume 
a simplistic surface density distribution, given by a power law, S(r) oc r~' p , with p = 1, as is 
likely for a disk with a constant viscosity parameter a in a steady state of viscous accretion 
(e.g, Hartmann et al. 1998). The dust inner radius is fixed at 0.5 AU. However, we do not 
accurately treat the inner rim in our simple dust model, and hence only report emission 
from beyond r = 0.7 AU where we believe our results to be accurate. We also calculate an 
additional disk model with an inner radius of 0.1 AU to test the sensitivity of line emission 
to this input parameter. We adopt a disk size of 200 AU, as is typical for observed disks 
(e.g., Andrews & Williams 2007). 

2.5. Thermal balance, chemistry and gas temperature 

The gas heating processes included are collisions with warmer dust grains, grain photo- 
electric heating by PAHs and small dust grains, X-rays, cosmic rays, photochemical reactions, 
FUV pumping of H2 and exothermic reactions. We include a modest chemical network (84 
species, ~ 600 reactions, see GH04) that involves the dominant chemical species of H, He, 
C, O, Ne, Mg, Fe, Si, Ar and S and solve for the gas chemistry (see Table |2] for adopted 
gas phase abundances of these elements). We have added Ne and Ar (since GH04), as they 
can dominate the cooling in the upper, X-ray heated layers of the disk (e.g., Glassgold et al. 
2007). We adopt similar procedures as Glassgold et al. for ionization of Ne (and Ar), where 
X-rays are assumed to doubly ionize the species, which then can recombine with electrons 
(or undergo charge exchange reactions with H) to form singly ionized species. We ignore 
ionizations to higher levels as the electron recombination rates to lower ionized states are 
typically high and most of the X-ray ionizations of the predominantly neutral species lead 
to the doubly ionized state (e.g., Maloney et al. 1996, Glassgold et al. 2007). H2 forms on 
grain surfaces, by reaction of H with H~, and by three-body reactions. The latter are espe- 
cially important for disks with a low cth- The chemistry in our models closely follows that of 
photodissociation regions or PDRs (Tielens & Hollenbach 1985, Kaufman et al. 1999), with 
some modifications for geometry and the high gas densities in disks (see GH04). Photodis- 
sociation of H 2 and CO by FUV and X-ray photons and their self-shielding is calculated as 
described in Tielens & Hollenbach (1985) and Maloney et al. (1996) and described in more 
detail in GH04. 



- 11 - 



Cooling processes include collisions with colder dust particles and line transitions of 
ions, atoms and molecules. Radiative gas cooling is usually dominated by transitions in the 
infrared of various chemical species, and we have used the approach outlined in Hollenbach 
et al. (1991) to estimate the background infrared radiation field (from dust re-radiation) 
and to calculate its effect on the level populations of the coolants. We refer the reader to 
GH04 for other details of the heating, cooling and chemistry. 

Gas line emission and photoevaporative flows typically originate in the upper layers 
of the disk, where the visual extinction to the star Ay < 1. The physics and chemistry 
included in our models therefore mainly focus on calculating the disk structure in the upper 
regions accurately. We do not consider some processes likely to be important in the very 
dense midplane regions of the disk, such as accretion heating (likely to be important in the 
midplane at r < few AU) for a standard MRI driven a-disk model and freezing of chemical 
species on cold dust grains (important in the midplane regions of the outer r > 5 AU disk, 
e.g., Aikawa et al. 1999, Markwick et al. 2002, Davis 2007). 



2.6. Disk structure 

We construct 1+1D models of the gas disk with the surface density distribution and 
radial extent of the disk specified. For simplicity, the gas and dust are assumed to be 
well- mixed, i.e., in the work presented in this paper, we assume a uniform gas-to-dust ratio 
throughout the disk and ignore the fact that settling may result in an increase of the gas/dust 
ratio with height. The gas pressure is set by vertical equilibrium between thermal pressure 
and gravity in the disk. Density and temperature are determined to be consistent with the 
prescribed surface density distribution, pressure, as well as chemical and thermal balance. 
Appendix A provides some details of the numerical scheme for obtaining a disk solution. The 
numerical models for the gas thermal, density and chemical structure are computationally 
intensive, and we therefore simplify our dust radiative transfer to a two-layer model (see 
Chiang & Goldreich 1997, D'Alessio et al. 1998, Dullemond, Dominik & Natta 2001, Rafikov 
& De Colle 2006). This simple approximation is fairly accurate in reproducing the disk 
structure as compared to full 2D calculations (e.g., Dullemond 2002). We use the correction 
factors prescribed by Dullemond, Dominik & Natta (2001) to account for the possibility 
of the disk becoming optically thin in the infrared, smoothly transit dust temperature from 
interior to the surface, and use the prescription by Rafikov & De Colle (2006) for determining 
the grazing angle to the disk surface. The vertical extent of the model disk is equal to the 
local radius, approximately 6-7 scale heights, where the visual extinction to the star is almost 
negligible. 
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The stellar EUV field will create an ionized HII region above the optical, FUV and X-ray 
heated neutral layer and we do not explicitly solve for the thermal structure of this ionized 
layer. Instead, we use the analytical formulae of Hollenbach et al. (1994) for the density of 
the ionized gas at the ionization front (IF) and assume a fixed temperature of 10 4 K for the 
ionized gas above this front. We ensure that our FUV and X-ray heating calculations extend 
to a height where the pressure is equal to that at the base of the ionized EUV-heated layer 
(i.e., the IF). The gas density in the EUV-heated layer is again determined from vertical 
pressure equilibrium, and ionization dominates gas chemistry in this layer. 

We take the analytic equation for the electron density n b (r) at the base of this layer (i.e., 
at the IF, below which lies the predominantly neutral disk) from Hollenbach et al. (1994) 



where 



n b {r) = n m (r/r g ) v (1) 

n fi<pEUv\ l/2 oc 1n5 ,1/2 ( M* V 3/2 _ 3 

nw = Ci (i^ ) = 2 5 x 10 [m; ) cm (2) 

"•. = ™(^)AU (3) 

and where C\ = 0.293 is a numerical correction factor, <Peuv — 1O 42 042 s _1 is the EUV 
luminosity of the star, a r = 2.57 x 10~ 13 cm 3 s _1 is the electron recombination coefficient for 
H + at 10 4 K, and r g is the characteristic gravitational radius where the sound speed in the 
HII region equals the escape speed from the disk/star system. The exponent p in Eq. [TJ is 
3/2 for r < r g and 5/2 for r > r g . 

The thermal pressure at the base of the EUV region (in units of nT) is then given by 
P b (r) ~2x 10 4 n b (r)cm~ 3 K, or 

A(r) = 5 x 10 9 ^ 2 (^) " 3/2 (£) " cm" 3 K (4) 

where p is 3/2 for r < r g and 5/2 for r > r g . 

In our disk model, we calculate the neutral gas density and temperature vertically, 
and insert the IF when the pressure becomes equal to Pb{r). At that height, we set the 
gas temperature to 10 4 K and the H nucleus density to n^(r). Because of the sharp rise in 
temperature across the IF, n and Ay drop abruptly across the IF. Above the base of this 
IF, the density drops as 

n{r) = n b (r)e- z2/2H{r)2 (5) 
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with the scale height H(r) given by H(r) = r g (r/r g ) 3 ^ 2 for r < r g and H(r) = r for r > r g 
(see Hollenbach et al. 1994). The vertical hydrogen nucleus column in the EUV layer is then 

/ M \ ~ 1//2 

N EUV ~ H(r)n b (r) = 3 x 1O 19 0^ 2 [t^J cm ~ 2 for r ^ ( 6 ) 

Note that this column is independent of r for r < r g , and, assuming oh = 2 x 10 _22 cm 2 /H, 
042 = 1, and M* = 1M , corresponds to a vertical Ay from the base to infinity of only 
~ 1.5 x 10~ 3 . The column from the base to the star is given by 

JV*~ rn h (r) = 3 x 1O 19 0^ 2 (^) ^ (£) ^cm" 2 (7) 

Therefore, at r = r g , the Ay from the base of the IF (or the base of the EUV layer) to the 
star is also ~ 1.5 x 10~ 3 for our standard dust model and 042 = 1, and M* = 1M Q . 

The different input parameters and assumptions of our fiducial disk model are summa- 
rized for convenience in Table [3j 



3. Results and discussion 

3.1. Disk structure and chemistry 

We describe the vertical structure of a fiducial model of a O.O3M disk around a 1M Q 
star. The stellar and disk parameters are given in Tables [TH3l The vertical density and 
temperature structure are shown as a function of height at a distance of 8 AU from the 
central star in Figure [TJ Gas temperature contours of the inner r < 10 AU region and 
the entire 200 AU disk are shown in Figures [2] and [3] respectively. Figures H] and [5] depict 
the heating and cooling and disk chemistry for the same vertical slice through the disk as 
Figure [D 

Figure [T] shows that the vertical gas and dust temperature structure allows three distinct 
layers to be defined in the disk, (i) the midplane layer, (ii) the intermediate layer and, (iii) 
the EUV-ionized surface layer. In the midplane layer (where Ay> 10 to the central star) 
the gas and dust temperatures are identical. At Ay~ 10 — 1, gas and dust temperatures 
deviate, to begin the intermediate layer. Figure [1] also shows an ionized surface layer in 
the uppermost regions of the disk at r = 8 AU, where gas temperatures rise to 10 4 K at 
z ~ 4 AU or Ay ~ 0.01 from the neutral side of the IF. We will focus mainly on the 
intermediate layers (between the midplane and the ionized layers of the disk) where gas 
and dust temperatures are unequal, because this is where infrared emission lines and FUV- 
generated photoevaporative flows originate, generally at Ay ~ 10 — 0.01 to the central star. 
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Figure [2] shows the gas temperature structure in the inner r < 10AU where most of the 
mid-infrared line emission originates, whereas Figure [3] shows the overall gas temperature 
structure for the entire 200 AU disk. 

Midplane layer. In the densest midplane regions of the disk, densities and opacities 
are high at all radii so that dust and gas are collisionally and radiatively coupled. Here, 
heating and cooling of gas is dominated by dust collisions as the high densities allow efficient 
collisional transfer of thermal energy. Dust and gas temperatures are almost equal and set 
by the nearly blackbody infrared radiation field. Other heating and cooling mechanisms are 
ineffective due to the large opacity (A v S> 1) along the line of sight to the star and to the 
disk surface (Figs. [T]and|l]). Due to the shielding of the disk interior, steady-state, simplified 
chemistry results in high abundances of molecular species such as H 2 , CO, O2 and SO2 in 
the midplane (Figure [5]). However, note that the dust midplane temperature falls below the 
freezing temperature of water ice for radii > few AU. Inclusion of ice formation on grains, 
presently ignored in these models, may result in different midplane chemistry and most of the 
elemental O may be incorporated in water ice, thereby depleting gas phase CO, O2 and SO2 
abundances, for example. Appendix B provides a justification for our assumption of no freeze- 
out in the intermediate and EUV-ionized surface layers and indicates conditions under which 
freezing of species on cold dust grains may become important. Detailed chemical models of 
disks which include ice formation have been developed by previous authors (e.g., Aikawa et 
al. 1997, Willacy et al. 1998, Aikawa & Herbst 1999, Willacy & Langer 2000, Markwick et al. 
2002, Willacy 2007). However, as discussed earlier, we ignore ice formation as the midplane 
regions are not important for either infrared line emission or for disk photoevaporation 
calculations, the two specific aims of our disk model calculations. 

Intermediate layer. At intermediate layers, that begin where Ay ~ 10 — 1 to the star, 
gas and dust begin to thermally decouple. With increasing vertical height and decreasing 
opacity to the central star, dust grains are directly heated by stellar optical photons to 
higher temperatures (the "superheated" layer of the two- layer dust models). As Ay to the 
star decreases, there is also increased penetration of stellar high energy photons that heat 
the gas. Grain photoelectric heating (due to PAHs and small dust grains) initiated by FUV 
photons, and penetration of ~ 1 keV X-ray photons that directly heat the gas, result in 
higher gas temperatures. Physics and chemistry in the intermediate layer are dominated 
by stellar photons. Gas-dust collisions become less important at the lower densities in the 
intermediate layer, and cooling is mainly by gas line radiation due to decreased opacity in 
the lines. Therefore, it is here that the mid-infrared gas lines largely originate, and can be 
observed in emission above the blackbody emission from the cooler midplane region. Figure H] 
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shows that at r = 8 AU, much of the X-ray and FUV heating of the gas is deposited at 
z ~ 1 - 2AU, where A v ~ 10 - 1. At this radius the [SI]25/xm, [OI]63//m, [NeII]12.8//m and 
[ArII]7yum are strong emission lines. The [OI]63/im line is one of the main coolants over 
much of the disk in the intermediate layer. In the inner regions of the disk (r < 20 AU) 
gas temperatures due to X-ray heating are high and [NeII]12.8/im and [ArII]7/im transitions 
become important in the cooling of the predominantly neutral gas (T > 1000K) that lies 
just below the EUV-heated layer. 

Figure [5] illustrates the chemistry of some of the important species at a disk radius of 
8 AU, and shows that photodissociation of molecules and photoionization of atoms takes 
place in the intermediate layer. The H 2 /H transition occurs where Ay is ~ 0.6 to the star, 
CO/C is at Ay ~ 1.6, and C/C + is at Ay ~ 0.4. The density of gas at the H/H 2 transition 
is ~ 10 7 cm -3 , with a locally incident stellar FUV field of Go ~ 10 5 , and FUV absorption is 
therefore dominated by gas opacity (e.g., Hollenbach & Tielens 1999) mainly due to neutral 
carbon and H 2 . At the warm temperatures characteristic of this radius, H 2 is additionally 
destroyed by the chemical reaction H 2 +0 leading to OH, which can be the dominant route 
at higher gas temperatures. The H/H 2 transition is therefore slightly deeper in the disk than 
the carbon ionization front. The conversion of C to CO takes place at higher Ay (deeper) 
in the disk, as the location of this transition is determined by where the dust opacity to 
FUV photons is ~ 1. We also note that in disks with reduced dust (low values of <th), lower 
formation rates of H 2 coupled with higher destruction rates due to the formation of OH for 
example, can lead to the formation of CO at higher z and lower Ay, even where the H 2 
abundances are low. This result was also noted by Glassgold et al. (2004) who found in 
their X-ray heated gas disk models that the C/CO transition often was at higher z compared 
to the H/H 2 transition. 

In the inner disk (r ~ 1 AU) the high densities and strong FUV field make heating due 
to H 2 formation and exothermic photoreactions important in the Ay~ 3 — 1 regions. Grain 
photoelectric heating is also significant at these Ay. X-rays dominate the heating both at 
lower (10 < Ay < 3) and at greater (1 < Ay < 0.03) heights. [OI] is the main coolant in the 
upper regions (Ay < 0.4) of the intermediate layer. Deeper in the disk, where the densities 
are high (n ~ 10 9 cm -3 ) and gas temperatures ~ 1000K, CO vibrational cooling dominates 
(1.6 < A v < 0.4). The CO/C and H 2 /H transitions are at A v ~ 1.2, with photodissociation 
by X-rays and FUV being important (e.g., Glassgold et al. 2004, Nomura & Millar 2005). 
At even higher densities and Ay, dust collisions and H 2 are the dominant coolants. 

In the outermost regions of the disk, densities are low (since £ oc r" 1 and scaleheights 
are larger), and a lack of collisional coupling can cause gas and dust temperature deviations 
to be significant. The effects of gas heating by FUV and X-rays are less important far 
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from the central source, and gas does not get warmer than the dust until higher values of 
z/r. Lower gas opacities result in efficient line cooling, gas temperatures often drop below 
dust temperatures, and dust collisions can be a heating agent in the intermediate layer. At 
r ~ 100 AU and at A v > 5, gas and dust temperatures are equal. Above these heights, 
the gas temperature first drops below that of the dust. Water cooling spikes at Ay~ 4, 
but is generally unimportant. In these intermediate layers, X-rays and dust collisions heat 
the gas up to heights where A v ~ 0.4, and CO rotational lines dominate gas cooling in the 
regions 5 < A v < 1. [01] takes over as the main coolant at A v < 1. Above A v ~ 0.4, 
grain photoelectric heating becomes significant, and causes gas temperature to rise above 
the dust temperature. The H 2 /H transition and CO/C/C + transitions occur high in the disk 
at Ay~ 0.03 and Ay ~ 0.5 respectively. 

The EUV layer. The uppermost layer in a disk is completely ionized by EUV photons 
from the central star. At r = 8 AU, the EUV layer is situated at a height z ~ 4.2 AU, where 
the density in the neutral gas just below the ionization front is n ~ 2 x 10 6 cm -3 (Fig. [TJ. 
The EUV layer at r = 1 AU is at z = 0.25 AU and n ~ 3 x 10 7 cm~ 3 , while at r = 100 AU, 
it is at z ~ 70 AU with the density n~4x 10 4 cm -3 in the neutral disk just below. We 
numerically calculate the ionization structure of the EUV layer and predict line luminosities 
produced in this layer in Hollenbach & Gorti (2008). 

3.2. Line emission 

Gas emission lines produced by dominant coolants can be important diagnostic probes 
of the disk conditions. Strong atomic and molecular cooling in the intermediate layers of 
disks where gas tends to be warmer than dust results in many mid-infrared emission lines 
originating principally from these layers. We discuss in this section the important transitions 
in the mid-infrared, far infrared and sub-millimeter wavelength regions and their detection 
with current and future observational facilities. These lines tend to be the most luminous 
lines arising from the ~ 1 — 200 AU disk, with the exception of some optical, near-IR and 
mid-IR lines produced in the very upper EUV-heated layer (Hollenbach & Gorti 2008) and 
some near-IR lines from thermally excited vibrational transitions of molecules in the inner 
(r < 5 AU) disk. Figure [6] shows the line spectrum from our fiducial disk of some strong 
emission lines in the infrared and sub-millimeter wavelength regions. Our disk models are 
computationally very intensive because a fine grid and many vertical iterations are necessary 
in order to resolve chemical transitions and for convergence of the disk structure solution. 
We therefore restrict ourselves to a limited parameter survey of four different models of disks 
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around a 1M Q star to study the dependence of the line strengths on the stellar radiation 
field and on the dust/gas ratio. Table H] lists the predicted line luminosities from our models 
for some of the strongest transitions. Model A is our standard run (Table [3]) which is also 
presented in Fig. El Model B has a dust opacity (<th) 100 times lower than Model A, and is 
considered to represent a more evolved disk that has had substantial grain growth or settling. 
Model C is a case with no X-rays and in Model D, the FUV luminosity of the star has been 
increased from Model A by a factor of 10. 

In the following, we discuss each of the strong lines and compare our results with model 
predictions by earlier authors when they exist. We also discuss our results in the context of 
the "Cores to Disks" (c2D) Spitzer Legacy Science Project data of gas line emission from 
young, optically thick disks (Lahuis et al. 2007). Spitzer IRS observations of disks by the 
c2D team were capable of detecting gas emission only from a few disks with very strong 
emission. The c2D survey of 76 circumstellar disks detected [Nell] emission from ~ 20% of 
sources, and [Fel]24/zm emission from ~ 9%. 8% of the sources had H 2 0-0 S(2) and/or S(3) 
line detections, and these were anti-correlated with the [Fel] detections. H2 S(l) emission 
was seen in only one source. Typical line luminosities range from 10~ 4 — 10~ 6 L for all 
the lines. We note that the many assumptions inherent to the models regarding stellar and 
disk properties make direct comparisons for each source possible only by individual detailed 
modeling. Here, in this first paper, we restrict ourselves to a more qualitative comparison of 
results with data, and model only one source, TW Hya, in detail. 



3.2.1. H 2 rotational lines 

The H 2 rotational lines arise from the inner disk, at r < 20 AU, and at heights where 
Ay to the star ranges from ~ 1 — 3 at the inner edge to ~ 0.3 at larger radii. The dominant 
H2 line in a region shifts as the temperature in the Ay~ 1 region drops with radius, from the 
9/im S(3) and 12/mi S(2) lines (at T~ 500 - 1000K) in the inner disk to the 17/im S(l) (at T 
~ 100 - 900K) and the 28/iin S(0) lines (at T~ 60 - 200K) from the 5 - 20 AU region. The 
S(3) and S(2) line emission peaks mainly at the inner edge (or rim, r < 1 AU)) where the 
gas temperatures are ~ 1000 K at these heights. S(l) emission is also significant here (30% 
of total luminosity), but most of the S(l) and all of the S(0) line emission is from r ~ 5 — 20 
AU. The S(l) line at 17/im is the strongest line for all our disk models, where we only report 
line emission from beyond 0.7 AU. Since most of the S(2) and S(3) emission arises from the 
inner edge of the disk in Model A, we calculated an additional case with a smaller inner disk 
radius. The S(2) and S(3) line strengths are nearly doubled when the gas disk is extended 
closer to the star (to 0.1AU), or if we consider direct illumination of the rim of a disk with 
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an inner hole (i.e., we include the 0.5 — 0.7 AU region in our line emission). In these cases, 
the S(2) line is comparable in strength to the S(l) line. 

Gas heating in the H 2 emitting region is mainly due to FUV-initiated grain photoelectric 
heating by PAHs, X-rays and FUV pumping of H 2 . In Model B, where the dust and PAH 
abundances are lowered by a factor of 100 from the standard Model A, gas heating due to 
the photoelectric effect is reduced and the decreased heating leads to lower gas temperatures. 
This is especially true in the r < 20 AU region. Here PAHs dominate the heating in the 
Ay~ 0.1 — 1 region and X-rays dominate at slightly higher Ay (Fig. Hj). Gas-dust collisional 
cooling is important in the Ay~ 0.1 — 1 region, and is lower relative to other coolants for 
Model B due to the lower <7h- The S(2) line arises from higher z (and lower Ay) where 
PAHs heat, and its strength decreases (at r < 5 AU) due to the lower gas temperature. On 
the other hand, S(l) emission increases here, which compensates for its decrease at greater 
radii (5 — 20AU), and the S(l) line strength remains relatively unaltered. The S(0) line is 
unchanged as it arises at greater Ay in the 5 — 20 AU disk where X-rays dominate heating and 
where the gas temperature is marginally affected by the reduced PAH abundance. Because 
of its lower excitation energy, the S(0) line is also less temperature sensitive than the S(l), 
S(2), and S(3) lines. The effect of neglecting X-rays (Model C) removes this heating agent, 
and the resulting lower temperatures decrease all H 2 line strengths in the disk. An increase 
in the FUV luminosity (Model D) by a factor of 10 from the standard model increases the 
strength of the lines as is to be expected. FUV heating in this case is important even at the 
higher Ay layers of the r < 20 AU disk region where the S(0) line originates. 

The rotational line luminosities of H 2 are ~ 10~ 5 L Q , well above the sensitivities of 
Spitzer and in favourable cases (luminous UV and X-ray stars that are nearby), ground- 
based facilities such as the TEXES spectrograph (Lacy et al. 2002) or MICHELLE (Glasse 
et al. 1997) on large telescopes. However, the high dust continuum at these wavelengths 
(~ 10~ 3 L Q for a resolving power of 700 provided by the IRS on Spitzer) indicates a very 
poor line-to-continuum ratio for Spitzer, and instruments with both high resolving power 
and high sensitivity are needed to see the lines. For very young disks around cTTs, with 
strong FUV heating due to perhaps a higher accretion rate (e.g., Gulbring et al. 1998), the 
H 2 lines are strong, and the S(2) line has been detected from several sources by the c2D 
team with the Spitzer Space Telescope. Our predicted H 2 S(2) line strength for a cTTs disk 
compares reasonably well with the typical observed value of 10~ 5 L Q (Lahuis et al. 2007). 
For computational tractability, we have restricted our disk models to radii larger than ~ 0.7 
AU. Extending our models to smaller inner disk radii results in thermal emission lines at 
higher excitation energies, such as the S(2) and S(3), being stronger. The strength of the 
S(0) and S(l) lines are comparable to the S(2) and S(3) lines in our models. However, S(0) 
and S(l) are undetected by the IRS, and we suggest that the decreasing line-to-continuum 
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ratio at longer wavelengths is perhaps the explanation for the non-detection. The c2D survey 
upper limits for the S(0) and S(l) lines are ~ 1O _5 L consistent with our model predictions. 
The EXES instrument on SOFIA has a sensitivity similar to that of the IRS, but a higher 
resolving power (R ~ 10 4 ) and thus is better suited to detecting lines in the presence of a 
strong continuum. It is capable of detecting lines with a luminosity of ~ 3 x 10~ 6 L for 
a source at the distance of Taurus (150 pc) and will probe gas in the 1 — 20 AU regions of 
nearby disks. However, limited sensitivity will make the detection of gas emission lines from 
distant disks challenging. 

Our H2 S(l) and S(2) emission line strengths are in agreement with the results of Nomura 
& Millar (2005), who obtain line luminosities of ~ 10~ 5 - 10~ 6 L for the H 2 S(l) and S(2) 
lines in their disk model with UV heating. Our results indicate that the FUV and X-ray 
luminosities of the central star are important in determining disk H 2 line emission, and that 
the dust properties (settling, coagulation) have a smaller influence. 

3.2.2. [Nell], [NelllJ and [Aril] 

The fine structure lines of [Nell], [Nelll], and [Aril] originate from high-temperature 
gas (T ~ 10 3 — 10 4 K) heated and ionized by X-rays and EUV photons from the central 
star. We discuss the emission from the EUV region in a separate paper (Hollenbach & Gorti 
2008) and here present only the calculations of [Nell] and [Aril] line emission from the X-ray 
heated neutral gas just below the completely ionized EUV region (Ay~ 0.02 — 2). X-rays 
ionize Ne and Ar in the neutral disk and the ionization fraction is ~ 10 — 20% in the emission 
regions (see also Glassgold et al. 2007). Radially, [Nell] and [Aril] emission from the X-ray 
heated gas is limited to < 20 AU, as first found by Glassgold et al. (2007). Beyond 20AU 
the gas temperature falls significantly below the excitation temperature, AE/k ~ HOOK, 
for [NeII]12.8yum. At the very surface of the neutral region, close to the EUV-heated region 
where the gas is fully ionized, and at the inner edge of the disk, [Nelll] is abundant. One 
important effect, noted by Glassgold et al. 2007, is the charge exchange reaction of Ne ++ with 
H, which lowers the Ne ++ abundance and thus the strength of [Nelll]. [Nelll] luminosities 
from the X-ray heated gas are therefore low and typical [NeIII]/[NeII] line ratios are < 0.1. 
[Aril] arises from higher temperatures (heights) than the [Nell] line. We find that [Nell] can 
contribute significantly to the total cooling of gas in these regions, at Ay< 0.1 within ~ 20 
AU. 

Predicted [Nell] and [Aril] line luminosities are moderately high ~ few 10 -6 L for all 
our disk models, and [Nell] has been observed around many young disks to date (Uchida 
et al. 2004, Geers et al. 2006, Lahuis et al. 2007, Pascucci et al. 2007, Espaillat et al. 
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2007) . Calculated and observed [Nelll] luminosities are lower than the [Nell] luminosity 
by a factor of ~10. However, X-ray produced [Nell] fluxes often fall short of the observed 
line fluxes, suggesting a contribution from the ionized EUV-heated layer as well. Since we 
ignore the EUV-heated layer, Model C with no X-rays shows no emission from these lines. 
However, dust properties and FUV luminosity in the presence of X-rays do affect the line 
strengths indirectly via their effect on the overall disk structure. Decreased dust opacity 
(or grain growth) results in a more "settled" disk with a smaller flaring angle, reducing the 
intercepted X-ray flux at the surface and thereby the [Nell] and [Aril] line emission. On the 
other hand, a higher FUV field (Model D) produces higher gas temperatures at the surface 
and greater disk flaring. The X-ray heated surface layer in this case subtends a larger solid 
angle. The increase in X-ray photon penetration in combination with the extra FUV heating 
of the gas results in higher line luminosities. 

Our Model A run is similar to the standard model of Meijerink et al. (2007), and the 
[Nell] luminosities are in good agreement. Meijerink et al (2007) do not include cooling by 
[Nell] which we find is important for the disk structure. Consequently, their gas temperatures 
in the surface [Nell] emitting layer are in general higher, which should produce higher [Nell] 
luminosities. However, they use a dust-determined density distribution that produces a 
surface subtending a smaller solid angle, and this decreases the amount of X-ray flux absorbed 
by the disk and lowers the [Nell] luminosity. Thus, the two effects partially cancel and their 
results are similar to our Model A. 

To further study the differences between our model and that of Meijerink et al. (2007), 
we compared our model results with their results for a disk with no accretion heating, using 
an almost identical model disk. In this run, we let the density distribution be determined by 
the dust temperature structure (as in their model) and used their assumed dust properties. 
We find that our [Nell] luminosities are, in fact, a factor of ~ 4 lower. This difference shows 
the importance of [Nell] and [Aril] cooling on the gas temperature. 

We note that the inclusion of [Nell] emission from the EUV-heated region would enhance 
the predicted total [Nell] line luminosities of all the models. For a typical EUV luminosity 
from a solar-type star of <pi = 4 x 10 41 s _1 , the calculated [Nell] line luminosity from the 
ionized HH-region above the disk is ~ 1O~ 6 L . The [NeIII]/[NeII] line ratio depends on 
the EUV spectrum and is higher for stars with a hard EUV spectrum (Hollenbach & Gorti 

2008) . 
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3.2.3. [Fel] and [Fell] 

[Fel]24/xm and [FeII]26/im line emission can be significant and detectable, even though 
they are seldom important disk coolants. [Fel] emission originates relatively deep at a vertical 
height given by Ay~ 3 — 10 to the star, and at r ~ 4 — 30 AU. At these vertical heights the 
dust and gas temperatures begin to decouple and the gas (T gas ~ 100 — 200K) gets warmer 
than the dust. As the Ay to the star decreases with height, Fe is rapidly ionized to form 
Fe + . [Fel] line luminosities are ~ 10~ 6 — 2 x 10~ 5 L in all our disk models. [Fel] arises 
from deeper layers and is strongly affected by the disk opacity (Model B). Decreased dust 
opacity in the disk increases the density at the Ay~ 5 layer whereas the gas, heated mainly 
by X-rays, is at the same temperature as in Model A, and the line strength increases. The 
[Fel] luminosity increases by a factor of 20, to 2 x 10~ 5 L , with a reduction in dust opacity 
by 100. With no X-rays incident on the disk (Model C), the [Fel] luminosity decreases due 
to the reduced heating. However, closer to the surface and the Fe/Fe + ionization front, there 
is some heating due to PAHs in this model which raises the gas temperature slightly above 
that of the dust, to produce [Fel] emission. An increase in the FUV field (Model D) results 
in higher gas temperatures, but on the other hand, also causes photoionization to Fe + . [Fell] 
increases with the FUV luminosity of the star and can be stronger than the [Fel] line (Model 
D). 

Our predicted [Fel] luminosities for all the cases considered here are at least a factor of 
5 lower than the most [Fel] luminous Spitzer IRS detections of the c2D team (Lahuis et al. 
2007), where the [Fel] luminosity when detected is ~ 10~ 6 — 10~ 4 L . From our analysis, we 
find that [Fel] emission is strongest when the disk dust opacity is very low, as in Model B. 
Further drastic reductions in dust opacity make the disk only marginally optically thick and 
the two-layer dust radiative transfer is no longer applicable. We note that our earlier models 
of optically thin disks predicted enhanced [Fel] emission from massive disks with low values 
of <7h (Gorti & Hollenbach 2004). These results are all suggestive of increased dust settling 
or grain growth enhancing the gas/dust ratio and decreasing <th in the r < 40 AU regions 
of disks where [Fel] is observed. We also note that our assumed gas phase abundance of Fe 
(10~ 7 ) may be an underestimate in disks. We will address these issues in future work on the 
modeling of the observed [Fel] emission in disks. 

3.2.4. [SI J 

The [SI]25.2/im line arises from the ~ 5 — 20 AU region in disks, at Ay~ 0.5 — 3. This 
line is strong in our model disk, although it is often not the dominant coolant in the region 
where it arises. In the sulfur emitting radial regions of the disk, heating is mainly due to 
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X-rays and FUV photons, and the main coolants are CO, [01], [SI] and dust collisions. We 
include detailed sulfur chemistry in our models, and assume a gas phase abundance typical 
of the diffuse interstellar medium, 2.8 x 10~ 5 (Savage & Sembach 1996). Sulfur is molecular 
(SO2) in the interior of the disk, is photodissociated to form S in the 1—3 region 

(column density ~ 10 22 cm -2 for our fiducial case), and is photoionized to form S + in the 
upper layers of the disk (Ay< 0.3). [SI] emission is reduced slightly for the reduced dust 
opacity disk (Model B) due to decreased heating by PAHs. As X-rays are an important 
heating mechanism in these regions, the no X-ray model (C) shows a reduction in the [SI] 
line luminosity as is to be expected. Increasing the FUV field of the central star increases 
the heating and raises the gas temperature to increase [SI] emission in Model D. 

Lahuis et al. (2007) do not detect [SI] emission from their disks, though they do detect the 
H 2 S(2) line, which we find to be comparable in strength. This could be because of the rising 
continuum from dust at the longer wavelength of the [SI] line, making line detection difficult. 
Sulfur may also deplete on grains in the emitting region, decreasing line emission. Watson 
et al. (2007) report the presence of strong [SI]25/im emission from two class I protostellar 
disks, but this emission may arise from shocks in the outflows around the protostars. We 
also note that our [SI] line luminosities are higher than predicted by Meijerink et al. (2007). 
They assume a much lower gas-phase sulfur abundance, and the trace amounts of carbon 
not locked in CO allows the formation of CS in their models in the regions where [SI] arises, 
thereby reducing the atomic S abundance. In our models the elemental S abundance far 
exceeds the available carbon abundance so that the formation of CS does not appreciably 
deplete atomic sulfur. Further, they do not include the effects of FUV radiation, which we 
find very important for determining the chemistry (ionization potential of S is 10.36 eV), 
and also for gas heating. Our gas temperature in this intermediate layer is typically higher 
than their models due to the inclusion of FUV. 



3.2.5. [Sill] 

[SiII]34.8/xm emission peaks slightly further out in the disk than many other mid-infrared 
emission lines which tend to arise mostly from the region inside of 25 AU. [Sill] arises from 
regions between 5 and 40 AU. It is relatively unaffected by changes in dust /gas ratio and 
X-rays, but increases significantly with FUV luminosity due to the combined effects of higher 
gas temperatures and a greater abundance of ionized silicon. This line falls in a relatively 
insensitive wavelength region of the IRS spectrometer on Spitzer, and also at a wavelength 
with high dust continuum emission, and is therefore difficult for Spitzer to detect. 
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3.2.6. [01] 

The [OI]63/im line is often the strongest cooling transition in disks, and because of its 
low excitation temperature (232 K) can arise from the surface regions of almost the entire 
disk. [01] emission originates in the atomic layer at 0.1 — 3.0 at all radii. The luminosity 
per logarithmic interval (dL/dhxr) rises nearly linearly out to ~ 20 AU, flattens beyond 20 
AU and then declines very slowly from ~ 120 AU. Because the line tends to be optically 
thick, the line luminosity is effectively determined by the gas temperature at the photosphere 
of the line, which typically occurs at Ay~ 1. In the regions where [OI] originates heating is 
due to both FUV photons and X-rays, and lowering the dust opacity (Model B) or the X-rays 
(Model C) only marginally decreases the line luminosity because the other mechanisms take 
over the heating. In Model D with higher FUV, there is more heating, the gas temperature 
is higher and [OI] emission is higher. Our calculated [OI] luminosities are high, ~ 10~ 4 L 
and this line is therefore a promising probe of the presence of atomic gas in the outer regions 
of disks. Our model results are also consistent with the results of Meijerink et al. (2007) 
who obtain an [OI]63yum line luminosity of ~ 5 x 10~ 5 for their standard model. The PACS 
spectrometer on Herschel is sensitive enough to detect OI emission from disks at 150 pc (e.g, 
Taurus, Chameleon and p Oph star-forming regions). PACS has relatively modest resolution 
and will not be able to resolve the line. The [OI] 145/im line is 10-100 times weaker than 
the 63/xm line, but may still be detected by PACS from disks around FUV-luminous sources. 
ISO detections of of the [OI] lines have been reported by Creech-Eakman et al. (2002) from 
the disks around AB Aur and other sources, and the observed luminosities ~ 10~ 4 L Q agree 
well with our model results. 



3.2.7. [CI] and [CI I] 

[CI] and [CII] are marginal coolants in the disk and increase in luminosity as the disk 
radius (and hence emitting area) increases. They originate where CO dissociates, and in 
general the C column density in the outer disk is low(~ 2 x 10 16 cm -2 at r ~ 100 AU), 
because the C is readily photoionized to C + . [CII] emission is stronger than [CI], and is high 
when the FUV luminosity is high (Model D) due to the higher gas temperatures and C + 
column densities in the outer disk. Our model results agree with the predictions of Kamp 
et al. (2006) and Meijerink et al. (2007) in their line luminosities, when we truncate the 
disks to match their assumed outer radii. The [CII]158/im transition, often a strong emission 
line along with the [OI]63/im line in photodissociation regions (e.g., Kaufman et al. 1999), 
is much weaker than [OI] in disks because of their higher densities and temperatures. The 
[CII] line may be too weak for detection by HIFI on Herschel. The sensitivity of the HIFI 
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instrument is about 10 times lower than that of PACS, and with the CII line luminosities 
calculated as being typically ~ 10~ 6 — 10~ 7 L , detecting [CII] emission from disks will be 
difficult, except for very nearby sources. 

3.2.8. CO rotational lines 

In the outer disk, CO is the main coolant from Ay ~ 5 where the gas temperature begins 
to deviate from the dust temperature, up to Ay ~ 0.5 where the CO is photodissociated. 
Often in this region, the gas temperature is less than the dust temperature, and thermal 
balance is achieved when the heating of gas by collisions with warmer dust equals the CO 
rotational line cooling. We also find that CO emission originates from heights where it is 
readily desorbed as ice from grains and is not significantly affected by freeze-out (Appendix 
B). In these intermediate layers, a reduction in dust opacity (Model B) marginally decreases 
the low J rotational line emission due to the lower gas temperatures. Low J rotational lines 
are relatively unaffected by X-ray heating, as X-rays are only important at depths greater 
than those characteristic of the surface layers producing the CO lines. An increase in the 
stellar FUV flux (Model D) can heat the gas in the outer disk (increasing line luminosities) 
but also increases the photodissociation rate of CO molecules, moving the CO surface to 
lower heights where gas temperatures are lower. We find that increasing the FUV field by a 
factor of 10 does not appreciably change the line luminosities. 

Ground-based CO observations with the JCMT, SMA and IRAM receivers are capable 
of detecting gas line luminosities of ~ 10~ 8 to 10~ 9 L from sources at ~ 100 pc. Our 
calculated low J CO line luminosities are much higher than these limits, and for nearby 
disks that are not significantly beam diluted, CO is a very good indicator of the presence of 
molecular gas in the outer (> 20 AU) regions. 

3.2.9. H 2 lines 

We estimate the strengths of H 2 lines using the derived density and temperature 
structure and solve for the level populations (lower 167 levels) after a disk solution for the 
temperature, density, and chemical abundance structure has been obtained using analytical 
approximations for the total cooling by H2O molecules (GH04). We note that cooling tran- 
sitions are often optically thick, and we use the escape probability formalism in calculating 
the line strengths as described in Hollenbach & McKee (1979). 

Water is a very important coolant in the disk in the lower parts of the intermediate layers 
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where the dust and gas temperatures begin to deviate (i.e., vertical heights where Ay~ 10 — 1) 
and where the gas temperatures are ~ 100 — 1000 K. Gas phase water abundances in the 
shielded interior of the disk are low, ~ 10 -8 -10 -6 , where H 2 formation is initiated by the 
weak ionization of H 2 to form which then reacts with H 2 to form Hjj". The Hjj" reacts 
with O to form OH + and subsequent reactions with H 2 ultimately leads to H3O" 1 ", which 
recombines with electrons to form H 2 0. At higher temperatures, typically > 300K, the main 
formation route is by neutral-neutral reactions of O with H 2 and OH with H 2 . The high 
temperatures are needed to overcome the large activation barriers of O and OH leading to 
H 2 0, and attained typically close to the Ay=l layer of the inner (< 40 AU) disk, where 
the hydrogen is still molecular. Water emission is strongest near the H 2 dissociation front, 
where gas temperatures and H 2 abundances are both high, and water is a strong coolant here 
(Fig. Hj). Line emission in the dense (n ~ 10 8 — 10 11 cm" 3 ) gas peaks as a function of z where 
H 2 and H 2 begin to photodissociate and gas temperature is rising. Although the total 
cooling is high (~ 3 x 10~ 5 L Q ), the emission is distributed over many different transitions 
at various wavelengths so that individual strong lines tend to be typically ~ 10~ 6 L Q . We 
note that in dense regions of disks, where the dust grain temperatures may fall below ~ 85K, 
water may freeze out on dust grains and may be depleted from the gas phase (see Appendix 
B). As we ignore freezing in our models, our predicted line luminosities may be upper limits 
for some of the longer wavelength and low-lying transitions of water. As shown in Appendix 
B, freeze-out of water may become important at large disk radii and at deeper Ay (where 
gas temperatures and photodesorption rates decline) where these longer wavelength emission 
lines originate. Table 5 lists some of the strongest infrared emission lines, many of which are 
detectable by HIFI and PACS on Herschel, by SOFIA (if atmospheric transmission permits), 
and a few mid-IR lines may be strong enough to be detected by Spitzer IRS. Line luminosities 
of transitions that may be affected by freeze-out are indicated as upper limits. The ortho- 
H 2 3i2 — 3o3 line at 274/xm, accessible by HIFI, is the strongest emission line and arises 
from the hot, dense gas at the inner edge of the disk (r < 1AU). H 2 line emission increases 
for lower dust opacities, because of higher gas densities at similar Ay, but only by a factor of 
~ 3 for a reduction in <th by 100 (see Table H]). In the disk with no X-rays, gas heating and 
destruction of H 2 are both lower, and the emission region shifts to higher z and slightly 
lower Ay where FUV heating is strong (see Fig. [4] for the Model A disk). The warmer gas at 
lower Ay but lower H 2 abundances act in opposition to keep the total water cooling nearly 
the same as the Model A disk. Increasing the FUV field (Model D) has similar effects of 
raising the gas temperature but simultaneously increasing H 2 photodissociation rates, and 
the total water cooling is not much affected. 
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4. Modeling the disk around TW Hya 

We apply our disk models to explain observed gas emission from the disk around the 
nearby, well-studied star, TW Hya. The properties of the central star and those of its dust 
disk are well determined, and therefore the number of free input parameters to our model 
are greatly reduced. Rotational CO emission line data from three different transitions have 
recently been obtained using the SMA by Qi et al. (2004, 2006) and we compare our model 
results with this uniform dataset which probes gas at different temperatures (and locations) 
in the disk. We also compare our results with the observed [Nell] emission by the Spitzer IRS 
(Uchida et al. 2004). This source is therefore well-suited for demonstrating the applicability 
of our models to determine the physical conditions and constrain the spatial location of the 
gas. 

TW Hya is the closest (51 pc, Mamajek 2005) known young T Tauri star (age ~ 4 — 10 
Myrs, Webb et al. 1999), and has a massive, face-on optically thick disk. TW Hya is known 
to be a strong X-ray source (Kastner et al. 1997). Muzerolle et al. (2000) estimate an 
ongoing accretion rate of 4 x lO~ lo M /yr and the star has an observed FUV spectrum (IUE, 
Valenti et al. 2000). TW Hya is also known to have a strong Lyman a emission component 
in its FUV flux (Herczeg et al. 2004). Lya excites the higher ro- vibrational levels of the 
ground electronic state of H 2 and populates higher electronic states. This is an important 
effect in determing the excited electronic and ro-vibrational emission of H 2 from the gas in 
the disk (Nomura Sz Millar 2005). For the largely thermally excited pure rotational emission 
discussed in this paper, the stellar Lya component is expected to be marginally important. 
For simplicity, we ignore the Lya emission from the star. The dust disk has been extensively 
modeled by Calvet et al. (2002) and extends from ~ 4 — 200 AU in size (Hughes et al. 
2007). Assuming a gas/dust mass ratio of 100, Calvet et al. determine the total disk mass 
(gas+dust) to be ~ O.O3M , if they take a maximum dust grain size of 1 cm. Table [6] lists 
the observational data and other input parameters used to model the gas emission. Our 
model disk is constructed by using the dust disk model of Calvet et al. (2002) and by adding 
a gas component with an interstellar gas to dust mass ratio (~ 100). We model the gas 
(and dust) components and calculate the density, temperature and chemistry as a function 
of spatial location in the disk. We compute the line luminosities of the CO and [Nell] lines 
and compare them with observations in Table [7] for two different disk radii. Model I has 
a radial extent of 174 AU, similar to that of the dust disk (e.g., Roberge et al. 2005) and 
Model II is a truncated gas disk with a radial extent of 120 AU. 

The [Nell] luminosity in the two models agrees extremely well with the observational 
[Nell] luminosity. In both models, the [Nell] line originates from the inner rim at 4 AU 
from X-ray heated gas exposed directly to the central star and so there is little dependence 
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on the outer disk radius. Though there appears to be good agreement between theory and 
observations, this result should be taken with caution. The inner rim is a region subject to 
flows, viscous accretion and photoevaporation (e.g., Chiang & Murray-Clay 2007, Alexander 
et al. 2006) and is not treated very accurately in our steady-state models. 

The CO lines in Model I are stronger than observed, a result also obtained by Qi et al. 
(2006) in their disk modeling. In our Model I, we overestimate the CO lines by a factor of 
~ 2. The CO 2-1 and 3-2 lines originate from the entire (174 AU) disk, whereas the CO 
6-5 line is emitted by slightly warmer gas confined to within ~ 100 AU from the star. The 
truncated disk Model II appears to better match observations, although we still overpredict 
the CO 6-5 line by a factor of ~ 1.6. 

We suggest that the gas disk around TW Hya may be truncated compared to the dust 
disk to ~ 120 AU, possibly due to photoevaporation. We calculate the disk photoevaporative 
mass fluxes (gE / dt) due to FUV radiation from TW Hya. TW Hya has a relatively high FUV 
luminosity ~ 4 x 10~ 3 L (IUE Atlas, Valenti et al. 2003). We use our disk temperature and 
density structure (Model I) and the analytical expressions of Adams et al. (2004) to estimate 
the photoevaporation timescales for the TW Hya disk. We find that the photoevaporation 
timescales (E / (dH / dt) where £ is the mass surface density of the disk) due to FUV radiation 
decreases with increasing disk radius and that the mass loss timescale at 120 AU is ~ 7 x 10 6 
years, similar to the age of the star/disk system. Therefore, it is plausible that during its 
lifetime, TW Hya has photoevaporated the disk outside of 120 AU, carrying with it all gas 
and small (< lOOyum) dust particles, but leaving behind the mm-sized dust grains seen in 
the sub-millimetre observations to extend to 174 AU. 

Table [7] lists other expected strong emission lines in our TW Hya disk model. We find 
that the [OI] line emission is very strong, well within the sensitivity limits of the PACS 
instrument on Herschel. PACS has a moderate resolving power (~ 1500) and hence will be 
unable to resolve the line, but is still capable of detecting the line (~ 10~ 5 L ) above the 
dust continuum at these wavelengths (~ 10~ 4 L Q at R~ 1500.). 

We end with a few caveats. Our models are steady-state and do not account for viscous 
spreading of the disk and photoevaporative flows which may introduce time-dependent or 
advection effects into the chemical and thermal structure. However, in the intermediate 
layer, the chemistry is largely driven by rapid UV photodissociation and chemical timescales 
tend to be shorter than dynamical timescales. Also note that our solutions may not be 
unique, and use of a different dust model (that is still consistent with the disk SED) may 
yield different results. We also note here that we do not require depletion of CO as has been 
inferred by earlier modelers (van Zadelhoff et al. 2001, Qi et al. 2004, 2006). We expect it 
to be unlikely that the CO is frozen on cold dust grains, as the grains in the intermediate 
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layer are sufficiently warm to thermally desorb CO ice mantles and the UV field of TW 
Hya is strong enough to keep the CO photodesorbed in the upper layers where the emission 
originates (also see Appendix B). 

5. Summary 

We present theoretical models of gas disks with optically thick dust around young stars. 
Our emphasis is on the physical and chemical processes relevant in determining the vertical 
density and temperature structure of the gas accurately in the surface regions, where Ay < 10 
to the central star, and where most infrared and sub-millimeter emission lines originate. 
We calculate expected luminosities of various diagnostic emission lines and find that many 
infrared lines are detectable by PACS and HIFI on Herschel and by EXES on SOFIA. Some 
lines are also detectable by the IRS on the Spitzer Space Telescope. Line luminosities are 
typically ~ 10~ 5 — 10~ 6 L . We find that decreased dust opacity in disks, a result of grain 
growth and/or settling processes during evolution, can increase [Fel]24/xm line strengths, 
but that other lines are marginally affected. X-rays heat the disk gas, and an absence of X- 
ray flux results in no [NeII]12.8/im and [ArII]7/xm emission from the predominantly neutral 
disk surface below the EUV-ionized layer. X-rays also penetrate deeper into the disk and 
are important heating agents in the Ay~ 1 — 3 regions, and affect the luminosities of the 
H 2 pure rotational lines, [Fel] and [SI] lines. The FUV luminosity of the central star is an 
important parameter, and in general, heats the gas to increase line emission. [Fell] and 
[Sill], in particular, increase with higher FUV fields. We find that the [OI]63yU.m line is a 
strong coolant and a good diagnostic probe of gas in disks. The high line luminosity in all 
our disk models, ~ 10~ 4 L , makes [01] a good signature of the presence of atomic gas and 
readily detectable by future facilities such as Herschel and SOFIA. 

We apply our models to the disk around TW Hya and can succesfully explain the 
observed CO and [Nell] emission from the disk. We suggest that the disk around TW Hya 
is being photoevaporated by its strong FUV radiation and that the gas disk is truncated 
(~ 120 AU) relative to the dust disk (~ 174 AU). 

6. Acknowledgements 

We thank the referee for useful comments that have improved this paper. We would 
like to thank Al Glassgold, Joan Najita and Rowin Meijerink for many helpful discussions 
during the course of this work, and for providing us with their results for model comparisons. 



-29- 



We would also like to thank Charlie Qi for providing us with his CO data on TW Hya, 
and Michael Kaufman for many useful discussions. We acknowledge financial support by 
research grants from the NASA Origins of the Solar System Program (SSO04-0043-0032), 
Astrophysics Theory Program (ATP04-0054-0083), and the NASA Astrobiology Institute. 

A. Numerical solutions of disk structure 

Our numerical scheme is as follows. We divide the disk radial extent into logarithmic 
intervals in r, with approximately 200 annuli for each model. We use the prescribed surface 
density power law to fix the surface density in each radial annulus. As our models are 1+1D 
solutions, the disk solutions outward of a given radial position do not affect the disk structure 
interior to that radius. However, column densities of gas and dust along the line of sight to 
the central star have to be calculated through the gas and dust in the inner radial zones, and 
we therefore proceed radially outward from the inner disk edge. The first initial solution of 
the disk vertical structure is obtained by solving the dust radiative transfer and setting the 
gas and dust temperatures to be equal. Starting from this initial guess for the density and 
temperature, we then solve for the chemistry and thermal balance in the disk and compute 
the gas temperature as a function of vertical height in the disk, z. We begin our solution at 
the midplane and advance to the surface, as this procedure ensures a more rapid convergence 
in the vertical density solution, which is constrained by the prescribed value for the surface 
density at that radius. 

Many vertical iterations are necessary for the density and temperature structure at a 
given radius in the disk to converge to a solution which gives the correct surface density 
and satisfies the condition for vertical hydrostatic equilibrium. The condition for vertical 
hydrostatic equilibrium specifies the gas pressure, which may not yield a unique solution for 
the gas density and temperature at a spatial location (r, z). We find all density-temperature 
combinations that result in the required gas pressure at each grid cell, and that satisfy the 
physical condition of decreasing density with height. For cells that yield multiple n — T 
solutions, we arbitrarily choose the higher density solution. We also note that some disk 
models may have convectively unstable regions, where we set the temperature constant with 
z until a stable solution is found at greater z. The density in the convectively unstable region 
is set by the pressure criterion. 

The gas temperature is determined by simultaneously solving chemistry and imposing 
thermal balance. Gas cooling is a strong function of the column densities of the different 
gas species to the surface of the disk. In calculating the escape probabilites, we assume two 
possible vertical directions perpendicular to the midplane of the disk as in some optically thin 
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transitions the escape probability from a height z through the disk midplane may be non- 
negligible. We use the gas escaping column densities calculated from chemical abundances 
in the previous chemical solution vertically through the disk. The vertical grid spacing is 
calculated in an adaptive manner and is not held constant. This is particularly important 
in order to resolve important transitions such as the H/H2 (which also changes the number 
density of gas) and C/CO transitions as one approaches the disk surface. We track important 
coolants such as CO and H 2 as abrupt gradients in their abundances can lead to sudden 
changes in gas temperature and erratic solutions. When the abundance of an important 
coolant begins to change with z, the grid spacing Az is made finer to resolve the transition 
and obtain a smooth solution. On an average ~ 200 vertical grid cells are required for 
each radial zone. Vertical convergence is obtained with an average of 10-20 iterations and 
convergence is more rapid in the outer disk. Our final disk solution is accurate to 2% in 
temperature and surface density. 



B. Justification of Assumption of No Freeze-out in Intermediate Layer. 

There are two desorption mechanisms that help clear grains of ice mantles in the inter- 
mediate layers of disks: (i) thermal desorption and (ii) FUV photodesorption. 

Thermal desorption. The grain temperature T gr at which ice mantles thermally desorb 
can be calculated by equating the rate of gas phase molecules hitting and sticking to a grain 
to the flux of thermally desorbing molecules from the ice surface (see Hollenbach et al. 2008). 

n sVs na 2 = u s e~ AE ^ ( (Bl) 



T 9 r = — A 7: Z .. m (B2) 



Here, n s is the gas density of species s, v s is its thermal speed, a is the grain radius, v s is the 
vibrational frequency of the adsorbed molecule, AE S is the binding energy of the adsorbed 
molecule, and A s is the area occupied by a single adsorbed molecule. From this we find 

AE s /k 
In [Av s /{A 8 n s v s )\ 

Using AE H20 /k = 4800K, AE co /k = 960K, v Ha0 = 7.8 x llPs -1 , z/ co = 4.2 x lO^s -1 , 
A H2 o — A C o — 10~ 15 cm 2 , n C o — 10 3 cm~ 3 , and n H2 o — 1 cm -3 (Hollenbach et al. 2008, 
densities of H2O and CO estimated from our models at r = 10 AU and Av=l), we find 
T gr (if 2 0) ~ 85K and T gr (CO) ~ 20K. In our models, the largest (coldest) grains are 85 K 
at 10 AU and 20 K at ~ 140 AU. Therefore thermal desorption prevents substantial freeze- 
out of CO in the intermediate layers of nearly our entire fiducial disk, and prevents H 2 
freeze-out only for r < 10 AU. 
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Photodesorption Hollenbach et al. (2008) discuss photodesorption of water ice in molec- 
ular clouds. Photodesorption will clear an ice mantle when the photodesorption rate from a 
grain covered with ice equals the sticking rate of water molecules from the gas. 

F FUV Yira 2 = n^ov^oiia 2 (B3) 

where Fpjjy is the FUV photon flux and Y ~ 3 x 10 -3 is the photodesorption yield from 
water ice (Westley et al. 1995, Andersson et al. 2006, Hollenbach et al. 2008). Fpuv is 
given approximately as 

Lfuv „-1.8Ay ( -da\ 
* FUV - T ~. o e l B4 J 

where L FUV = 10 31 L 31 erg s -1 is the FUV luminosity of the central star, and hvpuv — 1-6 X 
10~ n ergs is the typical energy of a photodesorbing FUV photon. From these equations, we 
find the critical radius r cr inside of which water ice mantles are cleared due to phoodesorption 

r cr = ( , , LPUVY ) V2 e"- A v (B5) 



Note that a z dependence enters because of the z dependence of uh 2 o, Ay, and vh 2 o- Taking 
i 2 o 



fiducial values of nn 2 o — 1 cm 3 an d vh 2 o = 3 x 10 4 cm s 1 we find 



i<2 

2000Lj( 2 ( ^^3 ) AU at A v = 1 (B6) 
= SOOLgf ( UH2 ° - ) ^ AU at A v = 3 (B7) 



A cm~ 3 . 

Clearly, r cr depends on nn 2 Oi which is a function of both z and r. However, it appears from 
the fiducial value of nn 2 o — 1 cm -3 taken from our model at r = 10 AU and at Ay= 1 that 
photodesorption should clear our grains of water ice to a depth of A v < 3. 

There is, however, an important caveat to these results. Water ice may form on grain 
surfaces by the reaction of adsorbed O atoms with adsorbed H atoms, which forms adsorbed 
OH, followed by the reaction of adsorbed OH with adsorbed H to form adsorbed H 2 (see 
Hollenbach et al. 2008). For this to occur, the adsorbed O atom must not thermally desorb 
before an H atom strikes a grain. This requirement translates to T gr < 36K or r > 50 AU 
for Ay~ 1. Beyond 50 AU one can use Eq. IB5I to estimate the critical radius inside of which 
photodesorption clears ice mantles, but one needs to use the larger of no or nn 2 o instead of 
n H 2 o- Since uq = 10 3 cm -3 at Ay= 1 and r ~ 100 AU, we find r cr ~ 60L^ 2 AU. Therefore, 
there is the possibility of water ice freeze-out for Ay> 1 and r > r cr ~ 60AU. This could 
deplete the gas phase oxygen not in CO and reduce somewhat the predicted luminosities 
of those lines whose greatest contributions come from these regions. We have tested the 
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sensitivity of our results to freeze-out by calculating the expected luminosity of all our lines 
from our disk models where we omitted any emission from the freeze-out regions, implied 
by Equations IB2I and IB51 We find that the longer wavelength transitions of water that arise 
mostly from low- lying states (at E/k < 100K) can be less luminous than in Table 5 by 
factors ranging from 2 — 10. The line luminosities for these transitions listed in Table 5 
should therefore be taken as upper limits, and we have marked them appropriately. 
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Table 1: Stellar Input Parameters For Fiducial 1M Q Case 

Mass 
Radius 
Temperature 
Luminosity 
Accretion rate 
Log L FUV (erg s^ 1 ) 
Log L x (erg s" 1 ) 
Log EUV (s^ 1 ) 



Table 2: Gas Phase Elemental Abundances 



Element 


Gas Phase Abundance 


H 


1.0 


He 


0.1 


C 


1.4 x 10~ 4 





3.2 x 10~ 4 


Ne 


1.2 x 10~ 4 


Mg 


1.1 x 10~ 6 


S 


2.8 x 10~ 5 


Si 


1.7 x 10~ 6 


Fe 


1.7 x 10~ 7 


Ar 


6.3 x 10~ 6 



Table 3: Fiducial Disk Model - Input Parameters 



Disk mass 


0.03 M* 


Surface density 


S(r) oc r" 1 


Inner disk radius 


0.5 AU 


Outer disk radius 


200 AU 


Gas/Dust Mass Ratio 


100 


Dust size distribution 


n(a) oc a -3 - 5 


Min. grain size a m i n 


50A 


Max. grain size a max 


20/im 




2 x 10~ 22 cm 2 /H 


PAH abundance/H 


8.4 x 10^ 8 



1Mb 
2.61 R 
4278 K 
2.34 L 

3.0 x 10" 8 M Q /yr 
31.7 
30.4 
41.6 
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Table 4: Predicted Line Luminosities (Log L/L ) 



Line 


Model A 


Model B 


Model C 


Model D 




(Fiducial) 


0.01 x o- H 


L x = 


10 x L FUV 


Aril 7fim 


-5.87 


-6.25 




-5.50 


H 2 S(2) 12/wi 


-4.82 


-4.97 


-5.10 


-4.40 


Nell 12.8/xm 


-5.54 


-5.83 




-5.28 


H 2 S(l) 17/xm 


-4.53 


-4.48 


-4.85 


-4.22 


Fel 24/xm 


-6.03 


-4.62 


-.6.41 


-5.94 


SI 25/im 


-5.39 


-5.55 


-5.84 


-4.9 


Fell 26/im 


-6.47 


-5.96 


-6.42 


-5.55 


H 2 S(0) 28/im 


-5.14 


-5.16 


-5.31 


-4.79 


H 2 (total cooling) 


-4.46 


-4.08 


-4.55 


-4.26 


Sill 35/im 


-5.87 


-5.84 


-6.01 


-4.94 


01 63yum 


-4.34 


-4.38 


-4.53 


-3.62 


CII 158/im 


-6.30 


-6.04 


-6.38 


-5.44 


CI 371/xm 


-6.60 


-6.42 


-7.14 


-6.27 


CO 2-1 


-6.98 


-7.19 


-6.96 


-6.89 


CO 3-2 


-6.55 


-6.79 


-6.42 


-6.33 


CO 6-5 


-5.92 


-5.84 


-5.96 


-5.72 
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Table 5: H 2 Emission Line Luminosities For Model Disk A 



Transition* 1 


A 1 / /TY1 1 




Transition 


\ ( 1 /m 1 
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— AnA 
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3^ 471 

OO. t: 1 X 
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71 ^38 

1 1 .OOO 
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U. ZjO 1 


P 


6 K 1 
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U. O 1 C/ 
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' 07 
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71 Q4fi 
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c; 
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-fi 38fi 
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75 37Q 


-o.oyjo 4. 
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u 24 
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36 21 2 


-6 381 

VJ . OO ± 
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3i 2 


78 740 

I (J . I t:W 


874 1 

"■Oil 4, 
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^41 


^74 


fi 128 

-U. liO 
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89 0^0 


(S 0^2 
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3Q 37Q 
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83 989 
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u 42 
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u 33 


3Q 3Q8 
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-P 934 

U . Z. Ot: 
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778 1 
-0. 1 1 4, 
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u 43 
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4fl 33fi 


-P 9 c iQ 
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8fi0 1 
J.ouu 4. 
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-fi 401 I 
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^40 


— 3ci 
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4Q 981 
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^23 


- 4.1 A 

^14 


1 ^9 40t 


Pi 1 fi^i 1 

-u.±u<j 4, 
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^41 


— 
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4Q 33fi 


-a. 1 i>o 4, 


p 


°30 


q 

— o 22 


1 3fi 4Q9 

1 0\J .^LZ/iLi 


-fi 4^4 1 





6 3 4 


5 2 3 


4Q 

ttc/.OOc/ 


fi 971 


p 


3l3 


— 2 02 


1 t94 

1 OO .OZiTT 


^ 728 1 
-0. ( zo 4, 


P 


533 


-4 22 


53.137 


-5.991 


p 


3 22 


-3ia 


156.193 


-6.150 | 


P 


4 3 i 


— 3 2 2 


56.324 


-5.728 i 





3o3 


-2l2 


174.620 


-5.947 | 


P 


4 2 2 


-3 13 


57.636 


-5.754 | 





2l2 


-loi 


179.523 


-6.070 | 





4 32 


— 3 2 i 


58.698 


-5.772 i 





221 


— 2 12 


180.482 


-5.934 | 





8is 




63.322 


-6.474 


p 


220 


-2 U 


243.969 


-6.223 | 


P 


808 


-7 17 


63.456 


-6.475 





3l2 


-2 21 


259.985 


-6.370 | 





625 


— 5u 


65.164 


-6.254 


p 


111 


-Ooo 


269.268 


-6.378 | 





330 


-2 21 


66.436 


-5.685 i 





3l2 


_ 3o3 


273.196 


-5.336 


p 


331 


— 2 20 


67.088 


-5.783 | 


p 


2o2 


-111 


303.447 


-6.372 | 


p 


524 


-4is 


71.066 


-5.991 | 


p 


2u 


— loo 


398.636 


-6.499 | 



a Thc symbols V and 'p' denote ortho-H^O and para-H^O transitions respectively. 
.[These lines may be affected by water ice formation on grains and are upper limits 
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Table 6: Model Parameters For Disk Around TW Hya 



Spectral Type 
Mass 

Stellar Temperature 
Stellar Radius 
X-ray Luminosity 
UV Spectrum 



K7 

0.77 M 
4000 K 

1 R (as in Calvet et al. 2002) 

2.3 x 10 30 erg s" 1 

IUE data from Valenti et al. 2003 



Dust Model 

Q"min 
Q"max 

Disk Inner Radius 
Disk Outer Radius 
Dust/Gas Mass Ratio 
Total Disk Mass 



From Calvet et al. 2002 

0.01 /im 

1 cm 

4 AU 

174 AU 

0.01 

0.03 M m 



Table 7: Comparison of Model Line Luminosities (in L Q ) and Observational Data. 



Model I Model II Obs. Data 

(R out =174 AU) (R o ^=120AU) 



CO 2-1 


1.8 x 10" 


-8 


8.4 x 10- 


-9 


a 7.9 x 10" 9 


CO 3-2 


5.1 x 10' 


-8 


2.4 x 10- 


-8 


a 2.4 x 10~ 8 


CO 6-5 


7.8 x 10' 


-8 


7.5 x 10- 


-8 


a 4.6 x 10" 8 


Nell 12.8/xm 


9.1 x 10" 


-6 


9.1 x 10' 


-6 


b 8.6 x 10- 6 


H 2 S(l) 17/im 


2.0 x 10" 


-7 


2.0 x 10- 


-7 




SI 25.2//m 


3.0 x 10- 


-6 


3.0 x 10- 


-6 




OI 63yum 


1.3 x 10- 


-5 


1.0 x 10- 


-5 




OI 145/xm 


3.3 x 10- 


-6 


2.0 x 10- 


-6 




CII 158/im 


9.9 x 10- 


-7 


6.7 x 10- 


-7 




CI 371/xm 


1.1 x 10- 


-7 


5.9 x 10- 


-8 





a Qi et al. 2006 
& Uchida et al. 2004 
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Fig. 1. — Vertical temperature and density profile at 8 AU. The gas temperature and density 
are shown by solid lines, and the mean dust temperature is depicted by the dashed line. The 
dotted vertical line shows the location of the height where the visual extinction Ay to the 
star is unity. 



-42- 



Log T 




50 100 150 

Radius (AU) 

Fig. 3. — Gas temperature contour plot of the fiducial model disk. 



012345 012345 
Z (AU) Z (AU) 

Fig. 4. — The dominant heating and cooling agents in the disk at a radius of 8 AU, as a 
function of disk height. PAH heating here refers to the sum of the heating by PAHs and by 
our grain size distribution, although PAHs tend to dominate this sum. 
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Fig. 5. — The chemistry of hydrogen, carbon, oxygen and sulfur-bearing species at a disk 
radius of 8 AU, with the abundances of dominant species as a function of disk height. Note 
that midplane (z < 1AU) chemistry may be subject to formation of ices on grains (refer 
text). 
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Fig. 6. — Line luminosities as a function of wavelength for the fiducial disk model extended 
inward to 0.1 AU (i.e., disk is from 0.1-200 AU). Only a few strong water and CO lines are 
shown in the figure. The [OI]63/im line, S(l), S(2) and S(3) pure rotational lines of H 2 , and 
certain H2O lines are some of the strongest emission lines from the disk. 



